{-# LANGUAGE DeriveDataTypeable   #-}
{-|
    Module      :  Numeric.ER.Real.Base
    Description :  arbitrary precision floating point numbers
    Copyright   :  (c) Michal Konecny
    License     :  BSD3

    Maintainer  :  mikkonecny@gmail.com
    Stability   :  experimental
    Portability :  portable

    This module defines an arbitrary precision floating point type and
    its operations.  It should be viewed more abstractly as an instance
    of 'B.ERRealBase' when used as interval endpoints.
-}
module Numeric.ER.Real.Base.Float
(
    ERFloat
)
where

import qualified Numeric.ER.BasicTypes.ExtendedInteger as EI
import Numeric.ER.BasicTypes.PlusMinus
import Numeric.ER.Misc
import Numeric.ER.BasicTypes
import qualified Numeric.ER.Real.Base as B

import Data.Ratio

import Data.Typeable
import Data.Generics.Basics
import Data.Binary
-- import BinaryDerive

--debugMsg = unsafePrint
debugMsg msg = id


{-|
A floating point number with a given but arbitrary precision represented by its 'Granularity'.

    * base: 2.
    
    * granularity specifies the bit-size of both the significand and the exponent  

    * special values: NaN, signed Infinity and signed Zero
    
    * no denormalised numbers
    
    * operations unify the granularity of their operands to the maximum 'Granularity'
     
    * Rounding is always towards +Infinity.  
      For field operations, the rounded result is as close as possible to the exact result.
-}
data ERFloat =
    ERFloatNaN -- any number / bottom
        { 
            apfltGran :: Granularity -- >0
        }
    | ERFloatInfty 
        { 
            apfltGran :: Granularity, -- >0
            apfltSign :: PlusMinus 
        }
    | ERFloatZero
        { 
            apfltGran :: Granularity, -- >0
            apfltSign :: PlusMinus 
        }
    | ERFloat
        {
            -- represents:
            -- sign * (1 + (mant/2^gran)) * (2 ^ exp)
            apfltGran :: Granularity, -- >0  granularity
            apfltSign :: PlusMinus,
            apfltMant :: Integer, -- 0 .. (2^gran - 1)
            apfltExp :: Integer -- -2^gran..2^gran
        }
    deriving (Typeable, Data)
    
zero = ERFloatZero 10 Plus
    
{- the following has been generated by BinaryDerive -}
instance Binary ERFloat where
  put (ERFloatNaN a) = putWord8 0 >> put a
  put (ERFloatInfty a b) = putWord8 1 >> put a >> put b
  put (ERFloatZero a b) = putWord8 2 >> put a >> put b
  put (ERFloat a b c d) = putWord8 3 >> put a >> put b >> put c >> put d
  get = do
    tag_ <- getWord8
    case tag_ of
      0 -> get >>= \a -> return (ERFloatNaN a)
      1 -> get >>= \a -> get >>= \b -> return (ERFloatInfty a b)
      2 -> get >>= \a -> get >>= \b -> return (ERFloatZero a b)
      3 -> get >>= \a -> get >>= \b -> get >>= \c -> get >>= \d -> return (ERFloat a b c d)
      _ -> fail "no parse"
{- the above has been generated by BinaryDerive -}
    
    
{-| normalisation

  * ensures that the components are within their regions
  
  * possibly turning the number into a zero or infinity
-}
normaliseERFloat :: ERFloat -> ERFloat
normaliseERFloat flt@(ERFloat gr s m e) 
    | m < 0 = 
        normaliseERFloat $ 
        ERFloat gr s (2*m + grmax) (e - 1)
    | m >= grmax =
        normaliseERFloat $ 
        ERFloat gr s ((m - grmax + (rndCorr s)) `div` 2) (e + 1)
    | e > grmax =
        case s of
            Plus -> ERFloatInfty gr Plus
            Minus -> minERFloat gr -- round upwards!
    | e < -grmax = 
        case s of
            Plus -> ulpERFloat gr -- round upwards!
            Minus -> ERFloatZero gr Minus
    | otherwise = flt
    where
    grmax = 2^gr
normaliseERFloat flt = flt

ulpERFloat gr =
    ERFloat gr Plus 0 (-2^gr)

minERFloat gr =
    ERFloat gr Minus (grmax - 1) grmax
    where
    grmax = 2^gr

maxERFloat gr =
    ERFloat gr Plus (grmax - 1) grmax
    where
    grmax = 2^gr

rndCorr Plus = 1
rndCorr Minus = 0

increaseERFloatExp e flt@(ERFloat gr s m eOld) =
    ERFloat gr s mNew e
    where
    mNew = 
        -grmax + ((m + grmax +  (rndCorr s) * (ediff - 1)) `div` ediff)
                           --  .^^^^^^^^^^^^^^^^^^^^^^^^^ round upwards
    grmax = 2^gr
    ediff = 2^(e - eOld) -- assuming e >= eOld
increaseERFloatExp _ flt = flt

decreaseERFloatExp e flt@(ERFloat gr s m eOld) =
    ERFloat gr s mNew e
    where
    mNew = 
        -grmax + ediff * (m + grmax)
    grmax = 2^gr
    ediff = 2^(eOld - e) -- assuming e <= eOld
decreaseERFloatExp _ flt = flt


apFloatExponent :: ERFloat -> EI.ExtendedInteger

apFloatExponent (ERFloatInfty _ _) = EI.PlusInfinity
apFloatExponent (ERFloatZero _ _) = EI.MinusInfinity
apFloatExponent (ERFloatNaN _) = EI.PlusInfinity -- includes infinity
apFloatExponent flt = EI.Finite $ apfltExp flt
        

setERFloatGranularity ::
    Granularity -> ERFloat -> ERFloat
setERFloatGranularity gr flt@(ERFloat oldGr s m e) 
    | gr > 0 =
        normaliseERFloat $ ERFloat gr s newM e
    | otherwise =
        flt
    where
    newM = 
        (m * (2^gr) 
          + ((rndCorr s)*(2^oldGr - 1))) -- round upwards!
        `div` (2^oldGr)
setERFloatGranularity gr f = f { apfltGran = gr } 
        
setERFloatMinGranularity ::
    Granularity -> ERFloat -> ERFloat
setERFloatMinGranularity gr flt
    | gr > oldGr = 
        setERFloatGranularity gr flt
    | otherwise = flt
    where
    oldGr = apfltGran flt
        
apfltGranularity = apfltGran

{-^ see the documentation of 'ERRealBase.getBaseMaxRounding' -}
apfltGetMaxRounding ::
    ERFloat -> ERFloat
apfltGetMaxRounding f@(ERFloatNaN _) = f
apfltGetMaxRounding f@(ERFloatInfty _ _) = f
apfltGetMaxRounding (ERFloatZero gr _) =
    ERFloat gr Plus 0 (-(2^gr))
apfltGetMaxRounding (ERFloat gr s m e) =
    ERFloat gr Plus 0 (max (e - (toInteger gr)) (-(2^gr)))

instance Show ERFloat where
    show = showERFloat 6 True False
    
    
showERFloat numDigits showGran showComponents flt =
    showERF flt
    where
    maybeGran gr
        | showGran = "{g=" ++ show gr ++ "}"
        | otherwise = ""
    showPM Plus = ""
    showPM Minus = "-"
    showERF (ERFloatNaN gr) = "NaN" ++ (maybeGran gr)  
    showERF (ERFloatZero gr pm) = showPM pm ++ "0.0" ++ (maybeGran gr)
    showERF (ERFloatInfty gr pm) = showPM pm ++ "oo" ++ (maybeGran gr)
    showERF f@(ERFloat gr s m e) =
        decimal ++ (maybeGran gr) ++ maybeComps
        where
        maybeComps
            | showComponents = "{val="++ show (s,m,e) ++ "}"
            | otherwise = ""
        decimal =
            showPM s
            ++ show digit1 ++ "." ++ (concat $ map show $ take numDigits digits)
            ++ (if dexp == 0 then "" else "e" ++ show dexp)
        (dexp, digit1 : digits) 
            | noLeadingZerosDexp == -1 =
                (0, 0 : noLeadingZeroDigits)
            | otherwise =
                (noLeadingZerosDexp, noLeadingZeroDigits)
        noLeadingZerosDexp = dexpBound - zerosCount
        (zerosCount, noLeadingZeroDigits) = 
            stripCountZeros 0 preDigits
        stripCountZeros prevZeros ds@(d : dsRest) 
            | d == 0 = stripCountZeros (prevZeros + 1) dsRest
            | otherwise = (prevZeros, ds)
        dexpBound -- upper bound of dexp: f/10^dexpBound < 1
            | e >= 0 = intLogUp 10 (2^e)
            | e < 0 = 2 - (intLogUp 10 (2^(-e)))
        preDigits =
            getDigits $ (abs $ setERFloatGranularity gran f) / (ten ^^ dexpBound)
        ten = setERFloatGranularity gran 10
        gran = 10 + (max (4 * numDigits) gr)
        getDigits ff =
            digit : digits
            where
            digit :: Integer
            digit = truncate ff
            digits =
                getDigits ((ff - (fromInteger digit)) * ten)
            

{-
    Beware: cannot use List.elem with ERFloat because of
    the intensional nature of Eq (eg ERFloatNaN /= ERFloatNaN).
-}
instance Eq ERFloat where
    (ERFloatNaN _) == _ = 
        False
        -- error "cannot compare NaN"
    _ == (ERFloatNaN _) = 
        False
        -- error "cannot compare NaN"
    (ERFloatZero _ _) == (ERFloatZero _ _) = True
    (ERFloatInfty _ pm1) == (ERFloatInfty _ pm2) = (pm1 == pm2)
    f1@(ERFloat gr1 s1 m1 e1) == f2@(ERFloat gr2 s2 m2 e2) 
        | gr1 < gr2 =
            (setERFloatGranularity gr2 f1) == f2
        | gr1 > gr2 =
            f1 == (setERFloatGranularity gr1 f2)
        | otherwise =
            s1 == s2 && m1 == m2 && e1 == e2
    _ == _ = False    

isERFloatNaN (ERFloatNaN _) = True
isERFloatNaN _ = False

instance Ord ERFloat where
    {- compare NaN -}
    compare a b@(ERFloatNaN _) =
        unsafePrint ("ERFloat: comparing NaN: " ++ show a ++ " vs. " ++ show b) EQ 
--        error $ "ERFloat: comparing NaN: " ++ show a ++ " vs. " ++ show b 
    compare a@(ERFloatNaN _) b = 
        unsafePrint ("ERFloat: comparing NaN: " ++ show a ++ " vs. " ++ show b) EQ 
--        error $ "ERFloat: comparing NaN: " ++ show a ++ " vs. " ++ show b 
    {- compare infty -}
    compare (ERFloatInfty gr1 pm1) (ERFloatInfty gr2 pm2) =
        compare pm1 pm2
    compare _ (ERFloatInfty _ Plus) = LT
    compare _ (ERFloatInfty _ Minus) = GT
    compare (ERFloatInfty _ Plus) _ = GT
    compare (ERFloatInfty _ Minus) _ = LT
    {- compare zero -}
    compare (ERFloatZero gr1 pm1) (ERFloatZero gr2 pm2) = EQ
    compare (ERFloatZero _ _) (ERFloat _ Plus _ _) = LT
    compare (ERFloatZero _ _) (ERFloat _ Minus _ _) = GT
    compare (ERFloat _ Minus _ _) (ERFloatZero _ _) = LT
    compare (ERFloat _ Plus _ _) (ERFloatZero _ _) = GT
    {- compare regular -}
    compare (ERFloat _ Minus _ _) (ERFloat _ Plus _ _) = LT
    compare (ERFloat _ Plus _ _) (ERFloat _ Minus _ _) = GT
    compare (ERFloat gr1 Plus m1 e1) (ERFloat gr2 _ m2 e2) 
        | e1 < e2 = LT
        | e1 > e2 = GT
        | gr1 == gr2 = compare m1 m2
        | otherwise = compare ((2^gr2)*m1) ((2^gr1)*m2)
    compare f1@(ERFloat _ Minus _ _) f2@(ERFloat _ _ _ _) =
        compare (-f2) (-f1)
        
instance Num ERFloat where
    fromInteger n
        | n == 0 = ERFloatZero (B.defaultGranularity zero) Plus
        | n < 0 =
            normaliseERFloat $ ERFloat gr Minus m e
        | otherwise = 
            normaliseERFloat $ ERFloat gr Plus m e
        where
        gr = fromInteger e
        e = max (toInteger (B.defaultGranularity zero)) $ (intLogUp 2 $ abs n) - 1
        m = (abs n) - 2^gr
    abs f@(ERFloatNaN _) = f
    abs f = f { apfltSign = Plus }
    signum f@(ERFloatNaN _) = f
    signum (ERFloatZero gr Plus) = setERFloatMinGranularity gr 1
    signum (ERFloatZero gr Minus) = setERFloatMinGranularity gr (-1)
    signum (ERFloatInfty gr Plus) = setERFloatMinGranularity gr 1
    signum (ERFloatInfty gr Minus) = setERFloatMinGranularity gr (-1)
    signum flt = 
        case apfltSign flt of { Plus -> 1; Minus -> -1 }
    negate (ERFloat gr s m e) = ERFloat gr (signNeg s) m e
    negate (ERFloatZero gr s) = ERFloatZero gr (signNeg s)
    negate (ERFloatInfty gr s) = ERFloatInfty gr (signNeg s)
    negate f@(ERFloatNaN _) = f
    {- addition -}
    f1 + f2 -- ensure equal granularity:
        | gr1 > gr2 = f1 + (setERFloatGranularity gr1 f2)
        | gr1 < gr2 = (setERFloatGranularity gr2 f1) + f2 
        where
        gr1 = apfltGran f1
        gr2 = apfltGran f2
    f@(ERFloatNaN _) + _ = f
    _ + f@(ERFloatNaN _) = f
    (ERFloatZero _ _) + f = f
    f + (ERFloatZero _ _) = f
    (ERFloatInfty gr Plus) + (ERFloatInfty _ Minus) =
        debugMsg ("ERFloat: infty - infty -> NaN\n") $ 
        ERFloatNaN gr
    (ERFloatInfty gr Minus) + (ERFloatInfty _ Plus) = 
        debugMsg ("ERFloat: -infty + infty -> NaN\n") $ 
        ERFloatNaN gr
    f@(ERFloatInfty gr s) + _ = f
    _ + f@(ERFloatInfty gr s) = f
    f1@(ERFloat gr s1 m1 e1) + f2@(ERFloat _ s2 m2 e2)
        -- equalise the exponents: 
        | e1 < e2 = f1 + (decreaseERFloatExp e1 f2)
        | e1 > e2 = (decreaseERFloatExp e2 f1) + f2
        -- ensure positive comes before negative: 
        | s1 == Minus && s2 == Plus = 
            f2 + f1
        -- opposite signs:
        | s1 == Plus && s2 == Minus && m1 == m2 =
            ERFloatZero gr Plus
        | s1 == Plus && s2 == Minus && m1 > m2 =
            normaliseERFloat $
            ERFloat gr s1 (m1 - m2 - 2^gr) e1
        | s1 == Plus && s2 == Minus && m1 < m2 =
            normaliseERFloat $
            ERFloat gr s2 (m2 - m1 - 2^gr) e1
        -- equal signs:
        | otherwise =
            normaliseERFloat $
            ERFloat gr s1 (m1 + m2 + 2^gr) e1
    {- multiplication -}
    -- ensure equal granularity:
    f1 * f2
        | gr1 > gr2 = f1 * (setERFloatGranularity gr1 f2)
        | gr1 < gr2 = (setERFloatGranularity gr2 f1) * f2 
        where
        gr1 = apfltGran f1
        gr2 = apfltGran f2
    -- NaN:
    f@(ERFloatNaN _) * _ = f
    _ * f@(ERFloatNaN _) = f
    -- Infty
    (ERFloatInfty gr _) * (ERFloatZero _ _) = 
        debugMsg ("ERFloat: infty * 0 -> NaN\n") $ 
        ERFloatNaN gr
    (ERFloatZero gr _) * (ERFloatInfty _ _) = 
        debugMsg ("ERFloat: 0 * infty -> NaN\n") $ 
        ERFloatNaN gr
    f * (ERFloatInfty gr s) = ERFloatInfty gr $ signMult s (apfltSign f)
    (ERFloatInfty gr s) * f = ERFloatInfty gr $ signMult s (apfltSign f)
    -- Zero
    (ERFloatZero gr s) * f = ERFloatZero gr $ signMult s (apfltSign f)
    f * (ERFloatZero gr s) = ERFloatZero gr $ signMult s (apfltSign f)
    -- regular
    f1@(ERFloat gr s1 m1 e1) * f2@(ERFloat _ s2 m2 e2) =
        normaliseERFloat $
        ERFloat gr s mNew (e1 + e2)
        where
        s = signMult s1 s2
        mNew = 
            m1 + m2 
            + ((m1 * m2 + (rndCorr s) * (2^gr - 1)) 
               `div` 2^gr)
    
instance Fractional ERFloat where
    fromRational rat = 
--        error "ERFloat: fromRational cannot be implemented reliably: use apfloatFromRational instead"
        (fromInteger $ numerator rat) 
        / (fromInteger $ denominator rat)
    f1 / f2 
        | gr1 > gr2 = f1 / (setERFloatGranularity gr1 f2)
        | gr1 < gr2 = (setERFloatGranularity gr2 f1) / f2
        where
        gr1 = apfltGran f1
        gr2 = apfltGran f2
    f@(ERFloatNaN _) / _ = f
    f1 / f2 =
        case apfltSign f1 of
            Plus -> f1 * (recip f2)
            Minus -> (- f1) * (recip (- f2)) -- rounding upwards!
    recip f@(ERFloatNaN _) = f
    recip (ERFloatZero gr s) = ERFloatInfty gr s
    recip (ERFloatInfty gr s) = ERFloatZero gr s
    recip (ERFloat gr s m e) =
        normaliseERFloat $
        ERFloat gr s mNew (-e)
        where
        mNew = 
            (- grmax * m 
             + (rndCorr s) * (grmax + m -1)) -- rounding upwards!
            `div`
            (grmax + m)
        grmax = 2^gr
        
        
apfloatFromRational ::
    Granularity -> Rational -> ERFloat
apfloatFromRational gran rat = 
    (setERFloatMinGranularity gran (fromInteger $ numerator rat)) 
        / (fromInteger $ denominator rat)
        
     
        
instance Real ERFloat where
    toRational (ERFloat gr s m e) =
        case s of
            Plus -> r
            Minus -> -r
        where
        r = (eOn2R) * (1 + mR/(grOn2R))
        mR = toRational m
        eOn2R = toRational $ 2 ^^ e
        grOn2R = toRational $ 2 ^ gr
    toRational (ERFloatZero _ _) = 0
    toRational f = 
        error $ "cannot covert " ++ show f ++ " to a rational" 
    
instance RealFrac ERFloat where
    properFraction (ERFloatNaN _) = 
        error "no integral part in ERFloatNaN"
    properFraction (ERFloatZero _ _) =
        (0, 0)
    properFraction (ERFloatInfty _ _) =
        error "no integral part in ERFloatInfty"
    properFraction f@(ERFloat gr s m e) 
        | e < 0 = (0, f)
        | s == Plus =
            (n, frac)
        | s == Minus =
            (-n, frac)
        where
        n = fromInteger $ 2^e + (m*(2^e) `div` 2^gr)
        frac = f - (fromInteger $ toInteger n)
        
    
instance B.ERRealBase ERFloat
    where
    typeName _ = "ERFloat (pure Haskell implementation)"
    defaultGranularity _ = 10
    getApproxBinaryLog = apFloatExponent
    getGranularity = apfltGran
    setMinGranularity = setERFloatMinGranularity
    setGranularity = setERFloatGranularity
    getMaxRounding = apfltGetMaxRounding
    isERNaN (ERFloatNaN _) = True
    isERNaN _ = False
    erNaN = ERFloatNaN (B.defaultGranularity zero)
    isPlusInfinity (ERFloatInfty _ Plus) = True
    isPlusInfinity _ = False
    plusInfinity = ERFloatInfty (B.defaultGranularity zero) Plus
    fromIntegerUp i = fromInteger i    
    fromDouble d
        | isNaN d = ERFloatNaN (B.defaultGranularity zero)
        | otherwise = (fromRational . toRational) d
    toDouble (ERFloatInfty _ s) = signToNum s * (1/0)
    toDouble (ERFloatNaN _) = 0/0
    toDouble flt =
        (fromInteger $ numerator rat) / (fromInteger $ denominator rat)
        where
        rat = toRational flt
    fromFloat f
        | isNaN f = ERFloatNaN (B.defaultGranularity zero)
        | otherwise = (fromRational . toRational) f
    toFloat (ERFloatInfty _ s) = signToNum s * (1/0) 
    toFloat (ERFloatNaN _) = 0/0
    toFloat flt =
        (fromInteger $ numerator rat) / (fromInteger $ denominator rat)
        where
        rat = toRational flt
    showDiGrCmp = showERFloat
    